#Study 1: Empathy and Party Sympathy
library(stargazer)
library(plot.matrix)
library(ggplot2)
library(broom)

load("Data/study_1.RData")
df <- data


#-----Eelco Harteveld like thermometer figure-----
list <- c("SP", "GL", "PvdD", "PvdA",  "D66", "VOLT", "CU", "CDA", "VVD","SGP", "JA21", "PVV", "FvD")
Sympathy_by_party_voted <- aggregate(cbind(SP, GL, PvdA, PvdD, VOLT, D66, CU, CDA, SGP, VVD, JA21, FvD, PVV) ~ Voted, data = data, FUN = mean, na.rm = TRUE)
Sympathy_by_party_voted <- data.frame(Sympathy_by_party_voted[,-1], row.names=Sympathy_by_party_voted[,1])
Sympathy_by_party_voted$mean_party_sympathy <- apply(Sympathy_by_party_voted[list], MARGIN =  1, FUN = mean)
Sympathy_by_party_voted <- Sympathy_by_party_voted %>% mutate(across(where(is.numeric), round, 3))
Sympathy_by_party_voted <- Sympathy_by_party_voted %>% mutate(across(where(is.numeric), ~ .x - mean_party_sympathy))
#Sympathy_by_party_voted <- Sympathy_by_party_voted %>% mutate(across(where(is.numeric), ~ .x *10))
Sympathy_Matrix <- data.matrix(Sympathy_by_party_voted[list,list], rownames.force = T)
RdBu<-rev(RColorBrewer::brewer.pal(n =11, name = "RdBu"))
plot(Sympathy_Matrix, col=RdBu, digits=2,
     main="In-Out Party Sympathy Evaluation",
     ylab="Party Voted for",
     xlab="Party Evaluated",
     cex=0.75,cex.lab=1.25, cex.axis=0.75, cex.main=2, cex.sub=2)

# This was manually saved as "In-Out Party Sympathy Evaluation.png" or Figure A1.

#----Descriptive Statistics Tables-----
library(vtable)
data_descriptives <- subset(data, select = c(empathic_concern, perspective_taking, fantasy, personal_distress,
                             absolute_inparty_sympathy, absolute_outparty_sympathy,range_party_sympathy, 
                             Distant_Inparty_favoritism, Distant_Outparty_favoritism,Spread_Affective_Polarization,
                             interest, left_right_placement, ideological_extremism, news_exposure,
                             GSL,LFT_CAT,OPL_Kort,SOCIALEKLASSE2016))

st(data_descriptives, title="Descriptive Statistics Tables", 
   labels=c("empathic concern", "perspective taking", "fantasy", "personal distress",
            "abs inparty sympathy", "abs outparty sympathy", "range party sympathy",
            "distant inparty favoritism", "distant outparty favoritism", "spread affective polarization",
            "political interest", "left-right ideology", "ideological extremism", "news exposure", "gender",
            "age", "education", "social class"), 
   summ=c('notNA(x)','mean(x)','sd(x)','min(x)','max(x)'),
   fit.page = '\\textwidth', out = 'latex', file="Tables and Figures/tableA1.tex")
# Note that some of the labels have been changed manually after saving this


#----Correlation tables-----
data_cor <- subset(data, select = c(empathic_concern, perspective_taking, fantasy, personal_distress,
                                  absolute_inparty_sympathy, absolute_outparty_sympathy,range_party_sympathy, 
                                  Distant_Inparty_favoritism, Distant_Outparty_favoritism,Spread_Affective_Polarization,
                                  interest, left_right_placement, ideological_extremism, news_exposure))
correlation.matrix <- cor(data_cor)
stargazer(correlation.matrix, title="Correlation Matrix")
# This was manually saved as Table A2.

#----Regression Models-----
#hypothesis 1 -> EC increases Inparty Liking
#hypothesis 2 -> EC increases Inparty Disliksing

#unweighted
in_party1 <- lm(absolute_inparty_sympathy ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                  left_right_placement + ideological_extremism + news_exposure + interest +
                  as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
        , data = df)
in_party2 <- lm(range_party_sympathy ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                  left_right_placement + ideological_extremism + news_exposure + interest +
                  as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                , data = df)
in_party3 <- lm(Distant_Inparty_favoritism ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                  left_right_placement + ideological_extremism + news_exposure + interest +
                  as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                , data = df)

out_party1 <- lm(absolute_outparty_sympathy ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                   left_right_placement + ideological_extremism + news_exposure + interest +
                   as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                , data = df)
out_party2 <- lm(Distant_Outparty_favoritism ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                   left_right_placement + ideological_extremism + news_exposure + interest +
                   as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                , data = df)

polarization1 <- lm(range_party_sympathy ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                      left_right_placement + ideological_extremism + news_exposure + interest +
                      as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                , data = df)
polarization2 <- lm(unweighted_spread_affective_polarization ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                      left_right_placement + ideological_extremism + news_exposure + interest +
                      as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                , data = df)
polarization3 <- lm(Spread_Affective_Polarization ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                      left_right_placement + ideological_extremism + news_exposure + interest +
                      as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                    , data = df)

stargazer(in_party1,out_party1,polarization1) # copy/pasted, Table A3
stargazer(in_party3,out_party2,polarization3) # copy/paste, Table A4

#weighted
in_party1w <- lm(absolute_inparty_sympathy ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                  left_right_placement + ideological_extremism + news_exposure + interest +
                  as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                , data = df, weights=weight)
in_party2w <- lm(range_party_sympathy ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                  left_right_placement + ideological_extremism + news_exposure + interest +
                  as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                , data = df, weights=weight)
in_party3w <- lm(Distant_Inparty_favoritism ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                  left_right_placement + ideological_extremism + news_exposure + interest +
                  as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                , data = df, weights=weight)

out_party1w <- lm(absolute_outparty_sympathy ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                   left_right_placement + ideological_extremism + news_exposure + interest +
                   as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                 , data = df, weights=weight)
out_party2w <- lm(Distant_Outparty_disfavoritism ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                   left_right_placement + ideological_extremism + news_exposure + interest +
                   as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                 , data = df, weights=weight)

polarization1w <- lm(range_party_sympathy ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                      left_right_placement + ideological_extremism + news_exposure + interest +
                      as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                    , data = df, weights=weight)
polarization2w <- lm(unweighted_spread_affective_polarization ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                      left_right_placement + ideological_extremism + news_exposure + interest +
                      as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                    , data = df, weights=weight)
polarization3w <- lm(Spread_Affective_Polarization ~  empathic_concern + perspective_taking + fantasy + personal_distress +
                      left_right_placement + ideological_extremism + news_exposure + interest +
                      as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
                    , data = df, weights=weight)

stargazer(out_party1w,in_party1w,polarization1w) # copy/pasted; Table A5 in paper
stargazer(out_party2w,in_party3w,polarization3w) # copy/pasted; Table A6 in paper

# Replicate SCK results study 1


#Make effect Size Figure ---------
betas <- c(tidy(in_party1)$estimate[2],tidy(out_party1)$estimate[2],
           tidy(in_party3)$estimate[2],tidy(out_party2)$estimate[2]*-1, # this DV was reverse-coded
           1.57, -1.41, # these numbers are taken directly from the SCK output
           tidy(in_party1)$estimate[3],tidy(out_party1)$estimate[3],
           tidy(in_party3)$estimate[3],tidy(out_party2)$estimate[3]*-1,
           -0.51, 1.17) # these numbers are taken directly from the SCK output
ses <- c(tidy(in_party1)$std.error[2],tidy(out_party1)$std.error[2],
         tidy(in_party3)$std.error[2],tidy(out_party2)$std.error[2]*-1, # this DV was reverse-coded
         0.59, 0.59, # these numbers are taken directly from the SCK output
         tidy(in_party1)$std.error[3],tidy(out_party1)$std.error[3],
         tidy(in_party3)$std.error[3],tidy(out_party2)$std.error[3]*-1,
         0.6, 0.67) # these numbers are taken directly from the SCK output

labels <- factor(rep(c("inparty","outparty"), 6),levels=c("outparty","inparty"))
labels2 <- rep(rep(c("Direct\nreplication","Conceptual\nreplication","SCK\nresults"), each=2),2)
labels2 <- factor(labels2, levels=c("Direct\nreplication", "Conceptual\nreplication","SCK\nresults"))

data <- data.frame(betas,ses,labels,labels2)
data$lo <- data$betas - 1.96*data$ses
data$hi <- data$betas + 1.96*data$ses
data$var <- c(rep("Empathic Concern",6),rep("Perspective-Taking",6))

fig.effects <- ggplot(data, aes(x=betas, y=labels)) +
  geom_point() +
  geom_errorbar(aes(xmin=lo, xmax=hi), data=data) +
  facet_grid(labels2~var) +
  geom_vline(xintercept=0, colour="red") +
  xlab("Effect sizes") + ylab("Liking") +
  theme(strip.text=element_text(size=12),
        axis.text=element_text(size=12))

ggsave(fig.effects, file="Tables and Figures/fig_effects.jpg", dpi=1000, width=6, height=4)


# Bayes scores not included in the paper -----
# library(rstanarm)
# library(bayestestR)
# library(BayesFactor)
# 
# m1_bay <- stan_glm(absolute_outparty_sympathy ~  empathic_concern + perspective_taking + fantasy + personal_distress +
#                      left_right_placement + ideological_extremism + news_exposure + interest +
#                      as.factor(GSL) + as.factor(LFT_CAT) + as.factor(OPL_Kort) + as.factor(SOCIALEKLASSE2016)
#                    , data = df)
# bayes1 <- bayesfactor_parameters(m1_bay, direction = 1)
# bayes_bcs5_social <- as_tibble(bayes1) %>% select(Parameter, BF) %>% rename("BCS70 Age 5" = BF)
# 
# m2_bay <- stan_glm(econ_scale ~ Restless + Fidgety + Fights + Worries + Solitary + Irritable + Miserable + Disobedient + Afraid_New + sex, data = df_final)
# bayes2 <- bayesfactor_parameters(m2_bay, direction = -1)
# bayes_bcs5_econ <- as_tibble(bayes2) %>% select(Parameter, BF) %>% rename("BCS70 Age 5" = BF)



